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■ Abstract. Some iterative techniques are defined to solve reversible inverse problems and a common 

formulation is explained. Numerical improvements are suggested and tests validate the methods. 

Resume. Nous definissons des techniques iteratives pour la resolution de problemes inverses reversibles 
, et en fournissons une formulation commune. Apres avoir suggere des ameliorations pour leur implementation, 

• ■ des experimentations sont presentees qui valident ces methodes. 



Introduction 



The classical Thermo Acoustic Tomography (TAT) problem formulates as follows. Consider an object /o 
contained in an open set Q of which emits an acoustic pressure wave at t = 0, considered as a Dirac pulse. 
This wave is modeled as a solution of: 



> 
m 
in 

O ■ f = in (0, T) X f2 



P(0) = /o, (1) 
dtP{0) = 0, 



where L is an operator modeling an acoustic wave phenomenon. Then this pressure wave is observed (e.g. 
^ thanks to piezoelectric sensors) and a set of observations is obtained from the solution p. That can be expressed 
• ^ ■ thanks to an observation operator C mapping a solution p to observations Cp. The inverse TAT problem 
rS ' consists in developing and studying methods to reconstruct /o from Cp and to define situations in which this 
^ reconstruction is possible. 

In the three past decades many techniques have been developed, offering effectual results (see works from 
authors of [1], ^ pTj, [13^, ^16^ among others). The new techniques we propose in section [1] rely on the following 
idea, influenced by [2]: if the system we consider is reversible in time, then the initial state to reconstruct can 
be seen, backward in time, as a state to reach, so that usual control and filtering techniques can be used to solve 
this inverse problem. For this purpose, we first used the Back and Forth Nudging algorithm (see [2 ) in |3]. 
With filtering techniques, as the Kalman filter defined in [7 and one of its reduced rank formulation, the SEEK 
filter (cf. [I2]), we deal here with possible improvements of this method. 

Of course, many assumptions are necessary to obtain a favorable observation situation: the way the wave 
propagates, depending on the media and on the kind of the wave, the final time and the number, size and 
position of the sensors recording the wave command the information contained in the data (see e.g. [S], [10] 
and the references therein). Moreover, even if the continuous problem is well set, numerical issues still put 
up some resistance, as considering noisy data, algorithmic complexity or apparition of spurious high-frequency 
oscillations during numerical implementations (the paper [TS^ surveys this point). 
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We deal here with this issue, as done in [13]. We introduce an artificial attenuation term in the numerical 
scheme that not only yields a regularization of the solution but corrects degenerated observation configurations 
for which filtering techniques are not helpful. 

1. Iterative stabilization and filtering for wave equations 

This section is devoted to main results about stabilization of the wave equation and Kalman-Bucy filter. 
Then we define iterative stabilizing methods. 

1.1. Observation and stabilization of the wave equation 

We assume that L = — A is the D'Alembert operator. If necessary, we can suppose that /o = and 
consider the stabilization problem for the initial value problem related to L instead of the inverse TAT problem. 
Similar considerations are also valid for linear variable speed (reversible) wave equations. 

Define p = ^, ^ , for any p e C^{{0,Ty,H^) fl CH(0, T); L^), and A = ^ on H = x L^, 

then D{A) = {H^ H H'^) x and equation ([1]) writes: 

P(o)^P.^(i°). 

The observations are defined in a Hilbert space U. It is convenient, when working with wave equations, to 
consider the time derivative of the observations, so that we assume that C G C{L'^^U) and use equally Cp' 
or Cp. In the practice of TAT, we only get observations from p and use its time derivative when needed. 

Concerning stabilizability and controllability of wave equations, we have the fundamental criterion: 

Definition 1.1. The observation inequality is satisfied if there exists T, M > such that: 

r \\Cp'\\ldt>M\\po\\l, (3) 

^0 

for all po G where p is the solution of ([T]) with initial data po- Scalar M is called observability constant. 

Indeed, in [TO], one finds the following result: 

Proposition 1.2. The three following propositions are equivalent: 
(i) The observation inequality ([3|) is satisfied. 

(ii) For every positive- definite self-adjoint operator T G L{U), the operator A — C^TC generates an expo- 
nentially stable Cq- semigroup on H. 
{Hi) The system {A^C^) is exactly controllable. 

Many geometrical interpretations to the observation inequality have been presented, mostly known as Geo- 
metric Optics Condition (GOC from [3j). In particular, these results explain this heuristic situation: when Cp' = 
l^p^ where u is an open subset of ft, then enough energy from po has to pass through u to get enough infor- 
mation to reconstruct po- It depends on many parameters such as the position of the sensors, the speed map 
of the wave equation, the final time T, etc. 

In this context, we introduce a first reconstruction method, the Kalman-Bucy filter. 

1.2. The Kalman-Bucy filter 

We recall the main results concerning the (continuous) Kalman-Bucy filter. It yields a way to approximate 
the real state that minimizes the error variance in the following situation (see [8j): 
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Assume that the true state solves the Hnear differential equation x*' = Mx* + v and the theoretical model 
is governed by x^ = Mx^. Given data y^, we denote the observation error by s, so as s = — Cx*. Errors e 
and V are null mean white Gaussian noise processes and their respective covariance matrices are Q and R. 

Definition 1.3. Given x^(0) and P^(0), the Kalman-Bucy filter consists in the two following differential 
equations, one to estimate the state x and one for the covariance matrix P, of a differential Riccati type: 

x' = Mx + PC^R-^(y^-Cx), 

P' = MP + PM^ + Q-PC^R-^CP^, 

and the Kalman gain K is given by K = PC^R~^. 

Let us explain some links between subsections II. II and II. 21 In the filters formulations, the feedback is realized 
thanks to an operator written PC^TC as feedbacks write C^TC and P = Idn in subsection 11.11 Thus one 
can consider here that P weights the feedback in comparison to the model. See ^17^ for some results similar to 
Proposition 11.21 about feedbacks PC^TC. 

We are studying different ways to define the stabilizing operator P, first with the nudging operator PC^TC = 
kC^C^ where /c > 0, that is in the framework of Proposition II. 21 then with filters. It leads to the following back 
and forth reconstruction algorithms. 

1.3. Iterative initial data reconstruction methods for inverse problems 

We get benefits from the reversibility of the wave equations, go back in time and use data again during a 
backward evolution, from t = T to t = 0, to deduce an approximation of /q. Such an idea has lead D. Auroux 
and J. Blum to define the Back and Forth Nudging algorithm in [2 , then D. Auroux and E. Cosme improved 
it with the Back and Forth SEEK (from private communication, see below for a description of the BF-SEEK). 
These techniques can be formulated for any sequences of positive operators (P/), (t/), (P|), (T|) as follows: 

Definition 1.4. Given a set of observations Cp^ and a rough estimate po, define Pq{0) = po. An iterative 
reconstruction method consists in iterating a back and forth process. The forward solution is given by: 

dupi = Api - PiC^Tidpi - p") in (0, T) x n, (4) 

with initial data p{(0) = pl_i{0) and p{'(0) = and the backward solution solves: 

dupi = Apl + PlC*nC{pl - p") in (0, T) X (5) 

with final data p\{T) = pl{T) and p\' {T) = pj^^T). The process is iterated for /c > 1. 

Concerning PC^TC, it is left implicit that time or space derivatives of (p{ — p^) and (p^ — p^) can be 
considered. The sign preceding P^C'^T^C{p^ — P^) changes for convenience of notation as one can notice that, 
when T contains a first-order time derivative, then the backward equation ([5]) writes forward: 

dupi=Apl-PiC*nC{pl-p% 

thanks to the variable substitution t ^ T — t. Finally, the correcting term still helps to stabilize the system. 

Note that unlike many usual methods, the model is considered here as a weak constraint, which can be useful 
since it may not be well known (e.g. refer to [B] about simulation of inhomogeneous acoustic speed model and 
relative issues in TAT). 

Finally, this kind of algorithms yields successive estimates {p\{^))k>Q of the initial object to reconstruct. As 
explained by Proposition ll.2l one knows that, under favorable observation conditions, both forward and backward 
equations are related to exponentially stable semigroups, that leads to the convergence of the algorithm. 
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2. Numerical experiments 

2.1. Discretization and numerical considerations 

Consider the 1-D domain Q = (-1/2, 1/2) and (0,T) = (0, 1) both uniformly gridded with steps 6x = 1/100 
and St = 1/200. We deal with the classical finite difference time domain theta-scheme (^-FDTD) to simulate 
the wave equation. Sensors are located periodically every Sdata grid points from the first one. 

Upper discussion about observability, stabilizability and controllability have their matching in this situation: 
a discrete observation condition occur, similar to ([3]), which is equivalent too to the discrete stabilizability of 
the system (we omit details, see [18 and the reference therein when ^ = 0). In order to compensate for possibly 
damaged observation conditions, due either to sensors configuration or to noise, we carry out the solution 
suggested in p^, adding an artificial viscous heating attenuation term in the 6>-FDTD, that leads to: 

Pn+l - 2pn ^ Pn-l _ .0 . . . Pn - Pn-1 p^^rr.^ ( Pn-1 - Pn- 
- A,Jp^_i,p^,p^+i) +£A,, ±PC TC ^ 

where A^^(pn-i,Pn,Pn+i) = OAsa:Pn-i + (1 - 2e)AsxPn^0A6a:Pn+u = 0.25, £ = or £ = ((5x)", a e (1,2), 
and ± stands to express both forward and backward implementations. 

The term ( Pn-i — Pn-i Pn — Pn ) allows us to consider derivatives of the correcting term in the feedbacks. 

2.2. From Kalman to SEEK filter 

Only main results about the Kalman filter are given and we describe then how to derive the SEEK filter from 
it (as in [15 ). These algorithms divide into two steps, a forecast one and an analysis one, in which are taken 
account the observations to correct the forecast, and two different kinds of parameters are considered, first the 
states (x^ and x^), then the relative error covariance matrices (P^ and P^) or their square root (S^ and S^). 
Definitions of the Kalman and SEEK filters are explained respectively on left and right columns below. 

f f f fl/2 

Definition 2.1. Given Xq, Pq and Sq = Pq , one iterates: 

Kalman Filter (KF) SEEK Filter 

Analysis step Analysis step 

K„ = [Pf-i + C^R-^C] C^R-\ K„ = Sf [Ir + {CSp^R-\CSi)] {CSi)^R-\ 

— ~ [C^n ~ Yn] ^ — ~ [C'x^ — y^] , 

Pn = - KnC] P^, S^ = S^ [ir + (CS^)^R"^ (CS^)] , 

Forecast step Forecast step 

Pj;+i = MP-MT + Q. Sj;+i = MS-. 

where Ik is the kxk identity matrix (in the discrete state space \ik = n and in the reduced rank space \ik = r) 
and Kn is the Kalman gain, minimizing the trace of the error covariance on xf , or the reduced rank Kalman 
gain in SEEK. 

Theoretically, one has P^ = S^S^^ + Q in SEEK, but since Q is not well known, we define P^ = (1+7)S^S^^, 
where 7 G (0, 1), which avoids an additional decomposition in SEEK. It is set similarly in KF. 

In KF, if the state space has a range of n, the forecast step necessitates 2n model evolution steps to get the 
forecast error covariance matrix from the analysis error covariance matrix, such that a less optimal gain may 
be considered to reduce calculation cost. This is the purpose of the SEEK to yield such a gain by considering a 
reduced rank gain more simple to obtain. Since the error covariance matrices are symmetric and positive-definite, 
Pham et al. [12] suggested to consider the following decomposition: if P is a symmetric positive-definite n x n 
matrix whose r < n larger eigenvalues are Ai > . . . > Ar > and their corresponding eigenvectors Vi, . . . , 14, 
the reduced decomposition of P is defined as the n x r matrix S = [y/XiVi . . . \^Vr]. Since SEEK consists in 
using the reduced decomposition of P-^, r model runs define the forecast error covariance matrix. 
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2.3. Results 

Four methods are compared: Time Reversal (TR), Back and Forth Nudging (BFN), Back and Forth SEEK 
(BF-SEEK) and Kalman Filter (KF). Let us set R = 0.3/, 7 = 0.01 for both BF-SEEK and KF, and the 
reduced rank r = 120, which is reduced when Rank [/ + (CSf )^R~-'^(CSf )] < r. An additive white Gaussian 
noise of level u is added to the data when said. The nudging gain is set to 0.9St for implementation limitations. 

Table [Hand Fig. [1] and [2] show RMS errors and some relative reconstructions obtained in various situations. 
The object to reconstruct is shown on the upper left part of figure [1] then, to the right, one sees TR and BFN 
reconstructions, followed by BF-SEEK and KF reconstructions. 



Settings 




TR 


BFN 


BF-SEEK 


KF 


Sdata 10 




5.4 


0.5 (100) 


0.6 (3) 


6 


^data = 10, 


TJ = 30% 


38.8 


16.3 (5) 


10.3 (2) 


15.7 


Sdata = 99 




8.3 


2.6 (100) 


6.3 (1) 


44.2 


Sdata = 99, 


TJ = 30% 


14.1 


10.8 (14) 


19.9 (1) 


46.4 


^data 99, 


u = 30%, a = 2 


9.8 


20.7 (8) 


12.2 (1) 


42 


Sdata > 99 




51.1 


20.7(7) 


83.3 (3) 
31.8 if a = 1.8 (3) 


95 



Table 1. RMS errors and number of iterations to converge (in brackets) for 10 {5 data = 10), 2 
{^data = 99) and 1 {Sdata > 99) sensor (s), noisy {u = 30%) or noiseless data, with or without 
the numerical attenuation Sx^. 




Figure 1. Object to reconstruct and reconstructions by TR, BFN, BF-SEEK and KF 
for 6data = 10 and u = 30%. 

As KF reacts quite well to noise addition, it shows much more sensitivity to the number of sensors and easily 
fails. BF-SEEK offers obvious improvements for both calculation cost and reconstruction error (Fig.[T]). When 2 
sensors are left, one can observe possible interesting effects of the attenuation term against noise with TR and 
BF-SEEK. Nevertheless, it may damage the reconstruction (with BFN) or have insignificant consequences (with 
TR). When only 1 sensor is left, only BFN and BF-SEEK are robust enough to yield a good approximation of 
the object, but BF-SEEK needs to be corrected with the attenuation to keep stable (Fig. [2j). 

3. Conclusions 

A common formulation for iterative stabilization of reversible evolution systems is given. It is used to define 
methods which solve some inverse problems for wave equations. Experiments show that the techniques we 
introduced may offer an alternative to usual inverse methods for TAT problem. In applications, knowledge 
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Figure 2. BF-SEEK reconstructions. Case Sdata = 99 and u = 30%: a = (a) and a = 2 (b). 
Case Sdata > 99: a = (c) and a = 1.8 (d). 

about the quality of the sensors, the background and the model approximations can be used by filters, but in 
this case they need to be precisely tuned. 

The use of an artificial attenuation is motivated by good results, and allows one to consider lossy medium in 
back and forth implementations when the physical loss does not exceed the numerical attenuation. 

When space dimension increases, we first face an excessive calculation cost (one BF-SEEK iteration with 
rank 60 equals almost one thousand BFN in computation time). So one would get interested in hybrid methods, 
e.g. getting first estimates with TR and BFN and then reconstructing the solution with filters. A solution is 
offered by Perfectly Matched Layers to reduce the space domain of calculation and for which a first order discrete 
scheme formulation is necessary for the filters. 
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